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Abstract 



We describe in detail the application of the recent non-Abelian Density 
Matrix Renormalization Group (DMRG) algorithm to the two dimensional 
t-J model. This extension of the DMRG algorithm allows us to keep the 
equivalent of twice as many basis states as the conventional DMRG algo- 
rithm for the same amount of computational effort, which permits a deeper 
understanding of the nature of the ground state. 
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Since the discovery of superconductivity in the rare-earth copper oxides there has been 
a growing interest in strongly correlated electronic systems. The t-J model proposed by 
Anderson (1987), and Zhang and Rice (1988) is an example of this interest. The theoretical 
predictions and implications of the model are possibly relevant and useful for a deeper 
understanding, particularly, the high temperature superconductors, and in generally the 
motion of "holes" in an antiferromagnet. 

It is conceptually a simple model and has been already widely studied in relation to the 
high temperature superconductors. However, it belongs to the class of systems which do 
not obey the condition of Perron- Frobenuis (Yosida 1980). This condition states that if the 
off-diagonal elements of a matrix are all non-positive and if the matrix is not in a block 
diagonal form then the ground state eigenvalue is non-degenerate. In the case of the t-J 
Hamiltonian the off-diagonal elements are not all non-positive. Thus the above theorem can 
not be applied, which implies that the phenomenon of ground state level crossing is present 
(Itoyama, McCoy and Perk 1990). As a direct consequence of this, the thermodynamic 
system will be unstable against phase separation. Indeed, as is known (Blatt and Weisskopf 
1979), in three-dimension, the system will collapse to arbitrarily high density. In two- 
dimensions we have the argument of Emery, Kivelson and Lin (1990) that the model is 
unstable against phase separation into hole-rich and no-hole phase. 

It is known that interchanging particles of different spin leads to a strong coupling be- 
tween the kinetic and spin degrees of freedom(Barnes 1989). Therefore, in different dimen- 
sions the model will represent different cases. Phase separation occurs also in one dimension, 
but for different values of the band filling and interaction strength than in the higher di- 
mension case. In two dimensions, it was argued that phase separation corresponds to stripe 
formation (Hellberg and Manousakis 1999). Stripe formation is one of the most controversial 
issues in the study of high temperature superconductors, where there is a phase separation 
of the holes which is limited to short range by, e.g.. Coulomb forces, and orders in striped 
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structures. Many experiments have found evidence for stripes Q and they arise in a number 
of different theories involving strong correlations (Zaanen and Gunnarsson 1989; Kivelson, 
Emery and Lin 1990; Zaanen, Horbach and van Saarloos 1996; White and Scalapino 1998). 

The first direct proof of the existence of a stripe ground state emerged from t-J model 
calculations (White and Scalapino 1998). The reason for this is that the t-J model inherits 
all the exchange hole correlations resulting from the antiparallel spin correlations in the 
f/ — s> oo limit of the Fermi sea. Hence, if stripes do exist, they will be more robust in the 
t-J model. Stripes have been studied by a number of numerical techniques in the t-J model, 
unfortunately resulting in conflicting conclusions. A major question is if stripes are part of 
the known phase separated regime of the t-J model (Hellberg and Manousakis 1999) or they 
represent a different ground state. This is the focus of our study. 

Using the recent developed non-Abelian DMRG (McCulloch and Gulacsi 2000) we study 
the ground state of the two dimensional t-J model in the vicinity of J/t = 0.35 and for dif- 
ferent size square lattices up to 24x6. We clearly find, for our choice of boundary conditions, 
a striped ground state. We have also noticed that in the phase separated regimes, the holes 
are attracted to the open boundary, while in the stripe regime, the holes are repelled by the 
boundary. This is evidence that that stripes are different from a usual phase transition phe- 
nomena, but also indicates that there are large finite size effects arising from the existence 
of the open boundary. 

The non-Abelian DMRG algorithm is numerically similar to the "Interaction Round a 
Face" (IRF) DMRG (Sierra and Nishino 1997). However, the decomposition of the Hamil- 
tonian into the blocks used by the DMRG algorithm is significantly different between the 
two algorithms. In IRF-DMRG, the vertex Hamiltonian is first transformed into an IRF 
model, and then a variant of DMRG is applied to the IRF Hamiltonian. In the non-Abelian 



For a recent review, see the proceedings of the Conference on Spectroscopies in Novel Supercon- 
ductors (Proceedings, 1998). 
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case, the standard DMRG algorithm is generahzed so that the block operators transform 
as irreducible representations of an arbitrary compact group (in this case SU(2)). This is a 
true generalization of the original DMRG algorithm (White 1992, 1993), in the sense that if 
the global symmetry group is chosen to be C/(l), as in the usual basis used in DMRG, the 
original DMRG algorithm is recovered exactly. In this paper, we fill in the details required 
to apply the non-Abelian algorithm to the 2 dimensional t — J model, in the basis given by 
particle number, N, and total spin, j, giving the global symmetry group C/(l) (8) SU{2). We 
then compare the results of this choice of basis with previous calculations using the usual 
basis, and comment on the conclusions that can be inferred from the DMRG calculations. 
The t — J Hamiltonian is 

H^-tY. i4aCja + H.c.) + J(^ Si ■ Sj - ^n,nj) (1) 

defined on the subspace of no double occupancy, and {i,j) means summation over nearest 
neighbour pairs only. The single site operators act on the usual three dimensional t-J basis of 
an empty site (hole), a single up spin and a single down spin. The difference for non-Abelian 
DMRG is that the spin up and spin down states form an SU (2) multiplet of dimension 2, thus 
reducing the single site basis to two states; a hole (transforming as the [0, 0] representation of 
U{1)®SU (2)) and a single spin (transforming as the [1, 1/2] representation of U{1)®SU (2)). 
All the single site operators are represented as 2 x 2 matrices acting on these states. The 
block basis states used in each step of the DMRG iteration are labelled \\nj{a)), which 
denotes the a'th state of n particles and total spin j. The matrix elements of the single 
site operators are simply the reduced matrix elements given by the Wigner-Eckart theorem 
(Biedenharn and Louck 1981), written here in our U{1) ® SU{2) basis, 

{n'fm'{a')\Ti,\njm{a)) = c2^,(ny(a')||T^||nj(«)) , (2) 

for the M'th component of an operator T"^ transforming as the [J] representation of SU(2). 
^mMm' denotes the Clebsch-Gordan coefficient. We have suppressed the trivial U{1) label 
on the T'^ operator. This theorem specifies the relationship between the reduced matrix 
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elements acting on the 2-state U{1) ® SU{2) basis, and the 3-state U{1) ® U{1) basis of 
particle number and 2;— component of spin f\ 

In matrix form, choosing basis vectors (1, 0) to be a hole, and (0, 1) to be a spin. The 
operators relevant to the t-J model are 





u 






, 





/ 


c ' ' = 


(0 
v 






,[0,1] = 


(0 











^[0,0] ^ 


(0 


0) 




p[o,o] ^ 







\ 






-1 





where we use square brackets to denote which representation of U{1) ® SU{2) the operators 
transform as. is the usual parity matrix used to enforce the correct commutation 

relations on the DIvIRG matrix operators. We can understand the unfamiliar matrix elements 
of the annihilation operator c^'^'^^^^ from the Wigner-Eckart theorem: when taking the 
hermitian conjugate of an operator T that transforms as some representation of SU (2), what 
we really mean is to find the operator such that (j'm'|T|/|jm) = (— l)"'~*'^(jm|T/^|jW)*. 
Due to the action of the Clebsch-Gordan coefficients in the Wigner-Eckart theorem, this is 
not the same as simply taking the Hermitian conjugate of the reduced matrix operator 
itself. Our choice is to make the full matrix elements of c'^^'^/^' identical with the usual 



^In a real calculation using the U{1) 03 ^(1) basis, the j index on the left hand side of Eq. 
cannot be used as a good quantum number since the usual DMRG density matrix in this basis 
does not preserve any non-Abelian global symmetries. 
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representation - the matrix elements of c''"'"^'^/^! then follow. The calculation of the matrix 
elements of the tensor product of two operators acting on different sites is given by the 
angular momentum theory as 
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SUi2) 



We have written this in a slightly convoluted way so as to make explicit the algebraic struc- 
ture. The coupling coefficients for f/(l) are so trivial that they are not usually represented 
in this form. Indeed, 



rii n2 n 
Ni N2 N 
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^n\,Ni+ni^ni,,N2+n2^n',n+N^n,ni+n2 
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(4) 



U{1) 



so the U{1) component of Eq. (H) reduces to exactly the coupling used in traditional DMRG. 
The point of departure from traditional DMRG is that the coupling coefficients for SU (2) 
are not so trivial: 



Ji 32 3 
ki k2 k 
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2j[ + lj2j!, + lJ2j + lV2kTl 
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31 32 3 
ki k2 k 



j'l 32 3' 
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where the {■ ■ ■} term is the Wigner 9j coefficient. 

We now write the Hamiltonian so that the operators transform as representations of the 
global symmetry group. For the t-J model, this is 



H 



-V2t Y: (c?'''/'1 ■ c^f'''^'^ + H.c.) - v^J E • 5f ■ ■ (6) 



By definition, the Hamiltonian itself transforms as the [0, 0] identity representation. We 
note that the interaction between two nearest neighbour sites in the non-Abehan 

representation requires summing 4 distinct terms, cjcj, qc], SiSj and niUj. In the basis, 
there are 8 terms: c|jC|j, c-i-jcj^-, cjjCjj, cj^icj^, S^S^, S^Sj, S~Sj' and riirij. Thus, although 
the matrix elements of the single site operators are more difficult to calculate using the 
non-Abelian formulation, there are correspondingly fewer matrix elements and operators 
required. 

We have applied this DMRG algorithm to the two dimensional t-J model by unroUing the 
two dimensional lattice into a one dimensional model with long range interactions, following 
the 'worm' approach (Liang and Pang 1994). We use periodic boundary conditions in the 
y direction and open boundary conditions in the (generally longer) x direction. Ideally, we 
would like to perform the calculations with periodic boundary conditions in both directions, 
but this would substantially increase the number of states required. Although the resulting 
one dimensional 'worm' model is reflection symmetric at the midpoint of the lattice, it is 
difficult to make use of this symmetry in DMRG calculations due to the non-uniform nature 
of the ground state. As a DMRG sweep progresses from one end of the system towards the 
centre point, the holes and spins tend to distribute themselves in a slightly asymmetric way 
between the left and right halves of the system, so that when the centre point is reached 
the left block basis is biased towards states that have too few holes, and the right block 
basis is biased towards states that have too many holes (or vice versa) . Enforcing reflection 
symmetry by using only one block plus its spatial reflection thus leads to a catastrophic 
reduction in the number of admissible superblock states, and a corresponding jump in the 
energy at that DMRG iteration. In all two dimensional models, constructing the initial state 
for the finite DMRG sweeps is more problematic than on one dimension. In one dimension, 
the so-called 'infinite' DMRG algorithm produces quite good results by itself, and usually 
provides a ground state that is quahtatively very similar to the converged finite DMRG 
state. In two dimensions however, the reduced translation symmetry of the unrolled one 
dimensional worm means that a genuine 'infinite' style algorithm would be rather difficult 
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^ While there are many possible ways to construct the initial blocks, we use the simple 
approach of constructing the initial blocks 'in place'; that is, starting from an initial 4 site 
system consisting of the 2 extreme sites from the left and right ends of the worm and adding 
two sites at a time, one from each end of the worm, working towards the centre of the system. 
This means that for most of the warm-up sweeps, there are no interaction terms between the 
left and right blocks. An alternative procedure is to rotate the system 90 degrees, so that the 
opposite ends of the worm are connected on the periodic boundary. However this introduces 
many more interactions between the left and right blocks throughout the calculation, which 
impacts on the accuracy. Another possibility is to introduce extra interactions for the warm- 
up sweep only, but we have not yet investigated this option. With no interactions between 
the two blocks, the eigenstates of the block density matrix coincide with the eigenstates of 
the block Hamiltonian. Thus there will only be a single non-zero density matrix eigenvalue 
(more, if the ground state is degenerate). The effect is that, until the first inter-block 
interaction appears, m — 1 of the block eigenstates are essentially random vectors. The 
initial state can be specified further by manipulating the target state as a function of system 
size. We have done this to obtain various initial condition: a state with all holes uniformly 
distributed, a phase separated state, and several random states. 

We have made calculations for various lattice sizes, keeping up to 1200 basis states per 
block. Table 1 shows a comparison of the ground state energy as a function of the number 
of basis states kept, using the U{1) Cg) U{1) and U{1) ® SU{2) basis, for a typical point 
in the 'striped' regime (White and Scalapino 1999; Xiang, Lou and Su 2001). The SU{2) 
symmetry provides a saving of a factor of two in the number of block states required. This 
is very significant as the computational complexity of the DMRG algorithm scales as at 
least 0{m^). However, even with 1200 states kept in the U{1) (8> SU{2) basis (equivalent to 
around 2500 states in the U{1) x f/(l) basis), the achieved energy is around 0.25% higher 



^Although there are interesting possibilities, see for example, Henelius (1999). 
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than the estimated true ground state energy. This compares very poorly with the accuracies 
generally achieved by DMRG for one dimensional models. 

In DMRG, the ground state wavefunction is iteratively improved, but only locally. This 
can lead to a situation where the DMRG converges self-consistently to an incorrect state, 
depending on the initial conditions and the details of the algorithm (Scalapino and White 
2000). In many cases, for a small number of states (but still relatively large compared with 
traditional DMRG studies) we have observed qualitatively different DMRG wavef unctions, 
depending on how we perform the warm-up sweep. We have also observed different converged 
wavefunctions even with the same initial condition, simply by varying the the rate at which 
the number of basis states per block is increased as the DMRG sweeps progress. For example, 
with the 16 X 6 system used in Table 1, using 500 basis states in the U (1) ® SU (2) basis, the 
ground state is most likely a two stripe configuration in agreement with (White and Scalapino 
1999). However, if we increase the number of retained states at a faster rate so that it takes 
fewer sweeps to reach the final total, we actually obtain a three stripe configuration. It is 
not until we increased the number of states to 800 that the three stripe configuration moves 
out of this local minima and formed the two stripe configuration. Simply doing more DMRG 
sweeps with 500 states is not effective. 

A possible way of dealing with this problem is to compare the energies of the competing 
DMRG ground states. The problem with this approach is that DMRG only provides a 
variational upper bound on the energy, and that the goodness of the variational energy (and 
therefore the truncation error associated with the DMRG state) can depend significantly on 
the nature of the ground state. Thus, this requires extrapolating the energies of both states 
to zero truncation error, showing that the difference in energy is statistically significant. This 
is a very time consuming procedure. In this paper, we only report results where we have a 
unique ground state, independent of the initial conditions (at least for the initial conditions 
we have used). However, the problem of multiple candidate ground states depending on 
the initial conditions deteriorates rather quickly as the system size increases. Indeed, one 
doesn't need to increase the system size too far before the DMRG fails to converge to a 
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believable ground state at all, at least in a reasonable number of sweeps. 

Fig. 1 shows the hole density along the x direction for a fixed number of holes, as the 
system size is increased. As the size grows, the stripes tend to move further apart while 
keeping approximately the same width, although it is difficult to make any real conclusions 
about the stripe width from this limited data. In particular, the width of the fluctuations 
in the real space density may be much larger than the correlation width of the stripe itself, 
if the stripe is delocalized. That depends on how much effect the boundary conditions have 
in pinning the stripes. Fig. 2 shows the effect of reducing the number of holes to 6, for the 
16 X 6 case. The ground state we observed is curious in that it breaks spatial reflection 
symmetry, but it provides evidence that it is energetically favourable to form a 'normal' 
stripe of hole density 4/6 and one 'proto' stripe, of hole density 2/6, rather than two stripes 
of equal hole density, or a single stripe. This suggests, as with the results from Fig. 1, that 
that the thermodynamic hole density per stripe is a constant (which depends on J/t). 

In Fig. 3, we attempt to flnd the optimal flUing per stripe, by increasing the system size 
and number of holes by 50%, to a 24 x 6 lattice with 12 holes. Since we consistently obtain 
three stripes, this hmits the hole density per unit length of the stripes to 0.5 < d < 1. 
This is the largest system that we could study, while still being reasonably certain that 
the obtained ground state is substantially independent of the initial conditions. Thus, for 
these parameters and boundary conditions, the two dimensional t-J model almost certainly 
has a striped ground state, and doping the system changes the density of stripes while the 
number of holes per stripe remains constant. It is difficult, however, to extrapolate these 
results to make definite conclusions about the nature of the ground state of the t-J model 
in the thermodynamic limit. Because of the half-periodic boundary conditions, the hole 
density is constant in the y direction. Thus any fluctuation in the hole density across the 
system, pinned by the open boundary in that direction, will appear as a vertical stripe in 
these calculations. Other possible groundstates of the thermodynamic t-J model, such as 
diagonal stripes or antiferromagnetic bubbles, are not permitted by construction. We also 
note that in the phase separated region, the holes are attracted to the open boundary of 

10 



the finite system. On the other hand, in the striped regime, the holes are repelled by the 
boundary. Thus the open boundary may well have a significant effect on the nature of the 
ground state that we observe, and especially, on the critical value of J/ 1 that separates stripe 
formation from phase separation. 

It is important to emphasize that stripe formation always implies the presence of an- 
tiphase boundaries, i.e., antiphase domain walls in the antiferromagnet. In all approaches, 
antiphase boundaries are always found in conjunction with stripes: from a simple mean-field 
calculations (Zaanen and Gunnarsson 1989) to the more sophisticated quantum numerical 
approaches (White and Scalapino 1998; Morais-Smith et al. 1998; Pryadko et al. 1999; Mar- 
tin et al. 2000). From our numerical data we cannot conclude that the antiphase boundaries 
are a consequence of stripes or vice-versa. What we are certain of is that the presence of 
the antiphase boundaries is a clear evidence of stripe existence. This favours our earlier 
observation that stripes are different from phase separation. 

An interesting explanation of the antiphase boundaries was suggested by Nagaev (1995) 
0. This brings us back to the general theory of metal-insulator transition, as formulated 
by Mott (1984), who pointed out that the number of free carriers should jump discontin- 
uously at the transition. Hence, there has to be a region of phase separation near the 
metal insulator transition. The existence of phase separation associated with doping away 
from an antiferromagnetic phase was recognized prior to the discovery of high temperature 
superconductors (Visscher 1974; Nagaev 1983). 

Nagaev refers to this phase as nanoscale phase separation (Nagaev 1995; Markiewicz 
1997), where phase separation is accompanied by charge separation, which in a perfect 
isotropic crystal, form an almost periodic structure (Nagaev 1995). Thus, a nanoscale phase 
separation is realized as a form of charge density wave. In this language the antiphase 
boundaries appear as the ferromagnetic droplets (Nagaev 1983, 1995): a hole can delocalize 



^For a review, see also Markiewicz (1997) and references cited therein. 
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over a finite domain by flipping the spins of the neighbouring Cu ion, forming a small 
ferromagnetic - ferron (Nagaev 1983, 1995) islands. It is interesting to note that if this is 
the origin of the antiphase boundaries then a second hole can lower its energy by localizing 
on the same ferron island, i.e., on the same antiphase boundary. This has the appearance 
of a real space pairing mechanism. 

Hence, looking at the stripes as being a nanoscale phase separation, the antiphase bound- 
aries (periodic ferrous) will appear as a consequence of the periodic hole structure (stripes) . 
At mean-field level this can be understood by recalling that the superposition (coexistence) 
of a charge and spin density-wave will always gives rise to ferromagnetism (Volkov, Kopaev 
and Rusianov 1973; Gulacsi and Gulacsi 1986). 

Work in Austraha was supported by the Austrahan Research Council. One of us (I. P. 
McCuUoch) acknowledges the hospitality of Los Alamos National Laboratory where much 
of the numerical work was performed. 
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TABLES 

TABLE I. Comparison of U{1) and SU{2) bases for the number of states versus ground state 
energy of a 16 x 6 t-J system with J = 0.35, t = 1, 8 holes and cyhndrical boundary conditions. 
The results using the U(l) basis are from reference (White and Scalapino 1999) We also include 
an estimate of the true energy, extrapolated to zero truncation error. 
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-52.65 ± 0.05 
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FIGURE CAPTIONS 



Figure 1. Hole density in the x direction for lattice sizes 16x6, 18x6 and 20 x 6, 8 holes, 
at J/t = 0.35. 

Figure 2. Hole density in the x direction for lattice size 16 x 6, 6 holes, J/t — 0.35. 
Figure 3. Hole density in the x direction for lattice size 24 x 6, 12 holes, J/t — 0.35. 
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Figure 1 , J/t=0.35, 8 holes 
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Figure 3, J/t=0.35, 24x6, 12 holes 
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